function out = getNonMarketUtility(outEstim)

outSim      = simProjection(outEstim.model,outEstim.setupFilter.shocks);

U0          = outEstim.model.params.U0;
U0d         = outEstim.model.params.U0d;
CHI         = outEstim.model.params.CHI;
B           = outEstim.model.params.B;
CHI0        = outEstim.model.params.CHI0;
PHIzero     = outEstim.model.params.PHIzero;
PHI         = outEstim.model.params.PHI;

posC        = find(strcmp(outEstim.model.labely,'$c_t$'));
posL        = find(strcmp(outEstim.model.labely,'$l_t$'));
posC_ba1    = find(strcmp(outEstim.model.labelx,'$c_{t-1}$'));
posD        = find(strcmp(outEstim.model.labelx,'$d_{t}$'));
posN        = find(strcmp(outEstim.model.labelx,'$n_{t}$'));
posZ        = find(strcmp(outEstim.model.labelx,'$\mu _{z,t}$'));

% the controls (includes the ss and risk-correction) and exp() to get the level
c_cu        = exp(outSim.y_cu(posC,:));
l_cu        = exp(outSim.y_cu(posL,:));

% the states - need to add the steady state
Css         = exp(outEstim.model.gSS(posC));
d_cu        = exp(outSim.x_cu(posD,:));
n_cu        = exp(outSim.x_cu(posN,:));
c_ba1       = exp(outSim.x_cu(posC_ba1,:)+outEstim.model.gSS(posC));
muz_cu      = exp(outSim.x_cu(posZ,:)+outEstim.model.hSS(posZ));

u_cu = U0 + d_cu.*(1/(1-CHI).*((c_cu-B.*c_ba1./muz_cu)./Css.^CHI0).^(1-CHI)+U0d ...
                + n_cu.*PHIzero./(1-1./PHI).*(1-l_cu).^(1-1./PHI));    

uNoConstant = d_cu.*(1./(1-CHI).*((c_cu-B.*c_ba1./muz_cu)./Css.^CHI0).^(1-CHI) ...
                + n_cu.*PHIzero./(1-1./PHI).*(1-l_cu).^(1-1./PHI));    

fracMarketUtility = uNoConstant./u_cu;
tmp = prctile(fracMarketUtility,[5 50 95.0]);
out.percentile5 = tmp(1,1);
out.percentile50 = tmp(1,2);
out.percentile95 = tmp(1,3);
out.mean         = mean(fracMarketUtility);
out.numSim       = size(u_cu,2);

end